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Abstract 

We calculate the bulk photovoltaic response of the ferroelectrics BaTiC>3 and PbTi03 from first 
principles by applying "shift current" theory to the electronic structure from density functional 
theory. The first principles results for BaTi03 reproduce eperimental photocurrent direction and 
magnitude as a function of light frequency, as well as the dependence of current on light polariza- 
tion, demonstrating that shift current is the dominant mechanism of the bulk photovoltaic effect 
in BaTiOs. Additionally, we analyze the relationship between response and material properties 
in detail. The photocurrent does not depend simply or strongly on the magnitude of material 
polarization, as has been previously assumed; instead, electronic states with delocalized, covalent 
bonding that is highly asymmetric along the current direction are required for strong shift current 
enhancements. The complexity of the response dependence on both external and material param- 
eters suggests applications not only in solar energy conversion, but photocatalysis and sensor and 
switch type devices as well. 
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7 INTRODUCTION 



The bulk photovoltaic effect - or photogalvanic effect - refers to the generation of intrin- 
sic photocurrents that can occur in single-phase materials lacking inversion symmetry jj— 14 ] . 
Ferroelectrics (materials that possess intrinsic, switchable polarization) exhibit this effect 
strongly, producing current in response to unpolarized, direct illumination. Traditionally, 
photovoltaic materials are heterogeneous, doped structures, relying on the electric field at 
ap-ti junction to separate photoexcited electrons and holes. By contrast, the bulk pho- 
tovoltaic effect can be observed even in pure homogeneous samples, as with BaTiC>3 js|. 
Recently, the effect has been demonstrated in the multiferroic BiFeC>3, with reported effi- 
ciencies as high as 10% 6|-|8|. Though ferroelectric photovoltaics are currently receiving a 
great deal of interest, the origins of their photovoltaic properties are considered unresolved. 
Attention has been focused on interface effects, crystal orientation, and the influence of 
grain boundaries and defects, while any bulk photovoltaic contributions have been largely 
ignored (q- 18] . Its mechanism has been proposed to be a combination of nonlinear optical 



processes, especially the phenomenon termed the "shift current" [19H22j . but this has not 
been firmly established, and the detailed dependence on material properties, especially in 
ferroelectrics, is not known. This has also hindered progress towards understanding other 
photovoltaic effects, as the bulk contributation could not be separated out. While shift 
current calculations have been performed for some non-ferroelectrics, no experimental com- 



26 parisons were performed 



23N25|. Here we present the first direct comparison of current 



computed from first principles with experimentally measured short-circuit photocurrent. 
Using the shift current theory, we successfully predict short circuit photocurrent direction, 
magnitude, and spectral features, demonstrating that shift current dominates the bulk pho- 
tovoltaic response. Additionally, we explore the relationship between material polarization 
and shift current response, making progress towards identifying the electronic structure 
properties that influence current strength. 
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We emphasize that nonlinear optical processes can give rise to a truly bulk effect 20, 



261 ] . The results demonstrate that the most important of these is the shift current, 



which arises from the second-order interaction with monochromatic light. The electrons are 
excited to coherent superpositions, which allows for net current flow due the asymmetry of 



the potential. Bulk polarization is not required; onl y in version symmetry must be 



bro 



Shift currents have been investigated experimentally 
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and computationally, though only for a few nonpolar materials 
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25j . Perturbative analysis 



treating the electromagnetic field classically yields the shift current expression 
J a = o'„ s E r Es 
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(^j) 2 E / dk - /Kk]) 6 (uv,(k) - uv(k) ± u 

n' ,n" 

90 n / n »(k, k) 



x (n'k| P r |n"k) (n"k| P s \n'k) 



dk n 



(1) 



34 where n',n" are band indices, k is the wavevector in the Brillouin zone, 0J n {k) is the energy 

35 of the nth band, so that a rs gives the current density response J to electromagnetic field E. 

36 Xrnj(k) denote the Berry connections for band n at k, and <f> n 'n" is the phase of the transition 

37 dipole between the bands n' and n" . It is worth noting that while the Berry connections 

38 introduce a gauge dependence, it is exactly canceled by the gauge dependence of d ' ' ' •" lVX} 

39 so that the overall expression is gauge invariant. 

We may view this expression as the product of two terms with physical meaning 



dk„ 



J dkr s (n',n",k;u)R g (n',n",k) 

n'.n" 



where 

r s (n',n",k;u) 
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J (/[n"k] - f[n'k}) (n'k| P r \n"k) (n"k\ P s \n'k) 

x 5 (u n » (k) — u n > (k) ± uj) (2) 

is the transition intensity, which is proportional to the imaginary part of the permittivity 
and describes the strength of the response for this transition, and 

d<fi n , n n(k,k) 



R q (n',n",k) 



dk n 



- \Xn"q{k) - Xn'g(k)] 



(3) 



is the shift vector, which gives the average distance traveled by the coherent carriers during 
their lifetimes. As an analytical tool, we compute and plot the quantity 



R q {uj) = / dk R q (n' ,n'' ,k)S (u n 'i(k) — u„'(k) ± u) . 



(4) 
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FIG. 1. The overall current susceptibility a, along with the imaginary component of the per- 
mittivity, e*, and shift vector integrated over the Brillouin zone R, are shown for [(a)] PbTi0 3 , [(b)] 
BaTiO-3. DFT-computed direct band gaps are marked with vertical lines. 



40 We note that R has units of length over frequency and is not physical, nor is it weighted 

41 by intensity; as such the R only provides qualitative information about the aggregate shift 

42 vector. For additional information, see the supplemental materials. 



43 METHODS 



44 Wavefunctions were generated using the Quantum Espresso and Abinit plane-wave DFT 

45 codes with the generalized gradient approximation exchange correlation functional. Norm- 



46 conserving, designed non-local pseudopotentials [38|, |39] were produced using the Opium 



47 package. Self-consistent calculations were performed on 8 x 8 x 8 k-point grids with energy 

48 cutoffs of 50 Ry; the resulting charge densities were used as input for non-self-consistent 

49 calculations on finer k-point grids as necessary. 



50 PbTi0 3 , BaTi0 3 



51 Lead titanate and barium titanate derive from the cubic perovskite structure and become 

52 tetragonal in the ferroelectric phase at room temperature with five-atom unit cells. Both 

53 exhibit strong, robust polarization; combined with their simplicity this makes them ideal 

54 candidates for investigating the structural influence on the shift current response. 



55 The calculations for 

56 perature geometries 
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bTiOs and BaTiOs were performed using experimental room tem- 



4 if ] . Shown in Fig. [T] are the shift current tensor elements, along 

57 with the imaginary component of the permittivity and Rz- Only current response in the 

58 direction of material polarization (Z) is shown. The two materials show broadly similar 

59 behavior, with the peak of response several eV above the band gap and well outside the 
eo visible spectrum, while the shift current at energies near the band gap is small. The shift 
ei current for both materials is stronger in response to incident light polarization parallel to 

62 the direction of ferroelectric distortion than when normal to it. 

63 The shift current depends weakly on the aggregate transition intensities and shift vectors. 

64 In PTO, at 5 eV there is a peak in both the current and intensity, yet the shift vector is 

65 relatively small. In fact, despite negative aggregate shift vector at many frequencies, the 

66 majority of the current response of PTO to xx polarized light is nonetheless positive. In 

67 BTO, the aggregate shift vector direction is largely positive with a negative shift current 
es response under zz polarized illumination. This indicates that contributions to response can 

69 vary significantly across the Brillouin zone, and suggests that strong correlations between 

70 large shift vector and high intensity response are possible but not guaranteed. The product 

71 of the aggregate transition intensity and shift vector does not determine even the direction 

72 of photocurrent; to find the shift current, it is vital to multiply the transition intensity by 

73 its associated shift vector, and then sum over bands and k-points. 



74 Experimental Comparison 



For bulk, single-crystal BaTiOs, experimental spectra are available for energies near the 
band gap ja, Q|- The total current in a bulk crystal for light incident normal to the current 
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direction can by computed from (J r „ s (oj) by 



J M = -^d E > ( 5 ) 

J q {u) = G r q r Iw (6) 

75 where a is the absorption coefficient, E is the electric field strength, which can be determined 

76 from the light intensity /, w is the width of the crystal surface exposed to illumination, and 

77 G rr is the Glass coefficient [3fl. This expression applies to samples of sufficient thickness to 

nn 

78 absorb all incident light. We obtained the light intensity and crystal dimensions from [5], I43| . 

79 which were ~ 0.4 mW/cm 2 and 0.1-0.2 cm, respectively. In Fig. [2} the experimental current 
so response from is compared to the response computed using shift current theory with 
si the parameters of the experiment (0.4 mW/cm 2 and 0.15 cm). Despite the uncertainty 

82 in experimental parameters, the agreement is striking, in both magnitude and spectrum 

83 profile, for both tensor elements. This includes the difference of sign between the majority 



84 of the transverse and longitudinal response, which is unusual [22|, as well as the small 

85 positive region of the longitudinal response near the band edge. For PbTiC>3, experimental 
se results suitable for quantitative comparison could not be located. However, we note that our 
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calculation for PbTiC>3 correctly predicts that the current direction is toward the positive 
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451 ] . as well as that it's 



88 material polarization for light frequencies near the band gap 

89 relatively insensitive to light polarization, in contrast to BaTiC>3. 

90 We emphasize that these calculations not only reproduce the magnitude of response, but 

91 its idiosyncratic features as well. Because this theory reproduces all the salient features 

92 found in the experiments, this comparison provides strong evidence that shift current is the 

93 correct description of the the bulk photovoltaic effect. 



94 POLARIZATION DEPENDENCE 



95 Presently unknown is the relationship of the bulk photovoltaic effect to the material po- 

96 larization. Identification of the bulk photovoltaic effect with shift current makes it clear 

97 that there is no direct, mechanistic dependence of response on material polarization, as 

98 is the case for many mechanisms to which photovoltaic effects in ferroelectrics have been 

99 attributed. However, shift current requires broken inversion symmetry, which here derives 

100 from the lattice distortion that produces ferroelectric polarization, suggesting that the re- 
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FIG. 2. For BaTi0 3 , the experimental current jj| and computed current (this work), for trans- 
verse (xxZ) and longitudinal (zzZ) electric field orientation, as a function of energy above their 
respective bandgaps. 



sponse may appear to depend on polarization in some fashion. However, Eq. ([T]) does not 
reveal a straightforward relationship between the magnitude of symmetry breaking, and the 
resulting shift current response. The presented data suggest that stronger polarization does 
not necessarily imply greater response; photocurrent densities in BaTiOs and PbTiOs are of 
similar magnitude, despite PbTiOs possessing more than double the material polarization 
of BaTi0 3 . 

To further investigate the connection of photovoltaic effect to material polarization, we 
studied a systematic family of structures based on PbTiOs. Starting with the cubic per- 
ovskite in the paraelectric structure, we rigidly displaced oxygen ions along a single Cartesian 
axis by amplitudes ranging from 0.01 to 0.09 lattice vectors, without otherwise altering the 
geometry. The spectra of shift current and aggregate shift vector are shown in Fig. [3] for 
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FIG. 3. The overall current susceptibility and aggregated shift vector R are shown for PbTiO-3 
with varying polarization. 



n3 several displacements. The results indicate a complex relationship between shift current and 
ii4 material polarization. As Fig3 shows, with soft mode amplitude 0.01, the shift current at 
us 3.2 eV above band gap is negative; with amplitude 0.07, the shift current reverses direction, 
lie resulting in a change of -200%. With amplitude 0.01, there is a negative peak at 3.8 eV; 
H7 with amplitude 0.07, the peak shifts to 4.2 eV and is four times the size, for an increase of 
us over 300%. However, in the intervening frequency range, the response is relatively small at 
n9 all displacements. 

120 Next, we turn our attention to the integrated shift vector R z (uS). The changes in shift 

121 vector are of special interest, since the symmetry constrains the overall shift current expres- 

122 sion via the shift vector. The integrated shift vector spectrum echoes the overall current 

123 response, but contains some distinct features. The increase in current from 4-5 eV does not 

124 appear to result from increased shift vector length, but from stronger coincidence of high 

125 transition intensity and large shift vectors. In fact, the overall shift vector changes little with 

126 displacement. However, from 7.5-8.5 eV, the integrated shift vector changes dramatically, 

127 suggesting that at some points in the Brillouin zone the oxygen displacement substantially 

128 alters the shift vector. Changes to the overall response are thus a combination of changes 

8 



129 both to shift vector and associated intensity. 

130 To understand these results, the electronic bands participating in transitions in these 

131 frequency ranges were examined directly. For the 4-5 eV range, examples of the transitions 

132 and associated Bloch states that dominate are shown in Fig. | ^a) at 0.01 and 0.09 lattice 

133 vector displacements. For this transition, the shift vector is 0.6 A at displacement of 0.01, 

134 and 1.0 A at 0.09. The valence state is largely composed of oxygen p-orbitals, while the 

135 conduction state is essentially a titanium d xy state. The states, like the shift vectors, are 

136 largely unchanged by the oxygen displacement. 

137 However, the transitions in the higher energy range are notably different. Shown in 

138 Fig. | ^b) are examples of the dominant transitions in the 7.5-8.5 eV range. The shift vector 

139 is large and positive (32.3 A) for 0.01 lattice vector displacement, and large but negative 
wo (-22.7 A) at 0.09 displacement. The participating valence state can be characterized as 
mi bonding between the Ti and O atoms collinear with polarization, while the conduction state 

142 features Ti-0 anti-bonding. These results point not to a simple dependence on material 

143 polarization, but to a dependence of shift current on the extent of localization of the initial 

144 and final states, which in turn depends on atomic displacement. Transitions between states 

145 that do not experience bonding interactions in the direction of ferroelectric polarization 
us manifest short shift vectors and insensitivity to oxygen displacement. 



us CONCLUSION 

149 The shift current response was calculated for ferroelectrics barium titanate and lead 

150 titanate. In the case of barium titanate, the shift current closely matches experiment, suc- 

151 cessfully predicting magnitude, sign, and spectral profile, including the notable dependence 

152 of current direction on polarization. For lead titanate, where reliable quantitative data is 

153 unavailable, the direction and its lack of dependence on polarization is reproduced. This 

154 strongly suggests that shift current is the dominant mechanism of the bulk photovoltaic 

155 effect in these materials. 

156 For the materials analyzed, the strongest responses are at frequencies well into the UV 

157 spectrum and outside the spectral range probed in most experiments. Consequently, the 

158 potential for large shift current response may not yet be fully realized, a conclusion supported 

159 by the very large bulk photovoltaic effect observed in response to X-rays |l9|]. Furthermore, 
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FIG. 4. (a) The non-bonding Bloch states of PbTi03 are involved in a transition that is insensitive 



to material polarization, with a shift vector length change from 0.6 A to 1.0 A as O sublattice 
displacement increases from 0.01 to 0.09, and |(b)| a transition from bonding to antibonding gives a 
shift vector that is highly sensitive to material polarization, with shift vector length change from 
32.4 A to -22.7 A for increasing O sublattice displacement. 
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160 the strength and direction of the photocurrent subtly depend on the electronic structure 
lei of the material, including covalent bonding interactions. This suggests that ferroelectric 

162 compounds can vary widely in response profile, and could potentially perform much better 

163 than previous results have indicated, encouraging efforts to design materials with large shift 

164 current response in the visible spectrum. 
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SUPPLEMENTARY NOTES 



Here we sketch the derivation of (jlj), as well as the discretized version suitable for com- 
putation 



Continuous Expression 

In the present work, the optical response is treated semiclassically. The quantum Liouville 
equation for the dipole Hamiltonian is expanded to second order. The perturbation and 
second-order density matrix term in the interaction picture are 



V(t) = — P • A coswt 
m 



h 2 



t r t' 



JO 



dt"dt' 



v(ty, v(ty,p(o) 



The shift current will arise from the (0; uj, —uS) terms, corresponding to successive interactions 
with the positive and negative frequency components of the incident electromagnetic field, 
creating a zero frequency, or non-oscillating, coherent state. The current density is 



J = - V [ dk (n'k'l p m \n"'k m ) (n"'k!"\ P 1 InV) 



\rnh) m 



dkS 



(n'k\ P ■ A \n"k) (n"k\ P ■ A \n"'k) (n'"k\ P \n'k) 

n'^n'" ,n" ,u 

x [{f[n"'k) - f[n"k}) 6 (w ft ,„(k) - w n „(k) ± co) 

1 



(f[n"k] - f[n'k}) 5 (cu n „(k) - u n ,(k) ± co)} 



oj n »,(k) - u n i(k) 



(1) 



The terms where n' = n'" cancel at k and — k due to time reversal symmetry, and are 
excluded. Thus, only coherent excitations can carry current. 
Using the identity 

(n'k'l P \ n k ) = (w n '"(k) — u n >(k)) (n'k'l X \n'"k'") , n" ^ n , k'" ^ k' (2) 



and the position operator for periodic systems first introduced in 42] 



(n'k'l X |nk) = -i8 nn ,V k 5(k - k') + 6(k' -k) I dx (n'k'\ zV k \nk) 

'wsc 



(3) 



we can obtain after some manipulation the expression 



Jq = cr r q s E r E s 



a r q s (u) = Tie 



x (n'k| P r \n"k) (n"k| P s |n'k) ( J^{k,k) _ [ Xn „ 9 ( k) _ Xn , g(k) ] 



(4) 



Discrete Expression 



Equation flU) cannot be calculated directly on a finite k-point grid. An appropriate expres- 
sion must be constructed that retains gauge invariance. Furthermore, it must accommodate 
the ambiguity in band identity at and near degeneracies. Adopting a similar strategy to [431 ] . 
we can reformulate the shift vector using 



Xq = l (n"k\^\n"k) 







x 1 1 dk' 



\xi{n"k\n"k') 



kg — kq 



= i lim — — In (n"k|n"(k + Ak q q)) - i\n (n"k\n"k) 

Ak q ^0 Ak q 1 

y = i lim — — In (n"k\n"(k + Ak n a)\ 



and 



(n'kl Pl-Hr- |n"k) ] = (n'kl P \n"k!) 



In (n'k\ P \n"k')) 



k> —k 



= (n'kl P \n"k) lim — (In (n'k\ P \n"(kAk g + q)) - In (n'kl P \n"k)) 
Ak^o Ak 

( „-ki P , » |„» k> ) = (rfk | P |„» k) Um ( h <"' k i * i""< k + Afc <i» 



(n'k|P|ra"k) 



Then we have 



P(n', n", k) = lim -——i I In 



(n"k\ P 


n 


'k> 


(n'k 


n'(k + Ak q q)) 




(n" 


k\P 


n'(k 


+ AA;,q)) 



+ ln 



(n"k| P 


n'k) 


(n"k 


n"(k + Ak q q)) 


(n"{k + Ak q q)\P 


n'k) 



2 



8 each term of which is manifestly gauge invariant. The expression can be generalized to the 

9 degenerate and near degenerate case by considering band sets and the appropriate operator 

10 sub-blocks. 



Glass Coefficient 



The solution to Maxwell's equations for and electromagnetic wave in a lossy medium can 
be written as E e lk ' r e luJt , where k is a complex wavevector, such that 

(k • k) = fiu 2 {iti + e r ) 

Solving for k(assuming propogation direction Z) 

k 2 z =fiu 2 (i€i + e r ) 



k z =± ^/Jiu 



y/e* + ef + e r / y/e* + e? - e r 



Also, 



/ oc \E\ 2 = E 2 e 2kiZ 



I = I Q e 2k * z 

which is Beer's law. Thus the absorption coefficient 



€7 — t,- 



a = 2ki = y/JIu^ 

can be determined from the permittivity, which from linear response theory is 
e(u) = e r (cu) + ieM = - f[m\) (m\ P \n) * V- 

mn 

+ i ^\\{f[n]-f[m\) (m\P\n) 

mn 

In a lossy medium, the intensity / depends on penetration depth z, and the total current 



10 - UJ m + u n 

2 

8iuj -u m + u n ) 



IS 



J q 


= «< / 


I rs dz = la 




J 




J q 


= ^7Io S 


[ e~ az dz 






Jo 


J q 


_.rs jrs 
°q J 

= w— 


[1 - e"^] 




a 



s e~ az dz 



where w is the width of the exposed surface, and I is the crystal depth. In the limit of 
al » 1 then 

^rs jrs 

J q = = wG r q s I r s (5) 



12 where G is the Glass coefficient l3i. 
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